Non equilibrium dynamics below the super-roughening transition. 
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The non equilibrium relaxational dynamics of the solid on solid model on a disordered substrate 
and the Sine Gordon model with random phase shifts is studied numerically. Close to the super- 
roughening temperature T g our results for the autocorrelations, spatial correlations and response 
function as well as for the fluctuation dissipation ratio (FDR) agree well with the prediction of 
a recent one loop RG calculation, whereas deep in the glassy low temperature phase substantial 
deviations occur. The change in the low temperature behavior of these quantities compared with 
the RG predictions is shown to be contained in a change of the functional temperature dependence 
of the dynamical exponent z(T), which relates the age t of the system with a length scale C(t): z(T) 
changes from a linear T-dependence close to T g to a 1/T-behavior far away from T g . By identifying 
spatial domains as connected patches of the exactly computable ground states of the system we 
demonstrate that the growing length scale C(t) is the characteristic size of thermally fluctuating 
clusters around "typical" long-lived configurations. 

PACS numbers: 



I. INTRODUCTION. 

Despite many efforts the understanding of non- 
equilibrium dynamics of disordered and glassy systems 
in finite dimensions remains a challenging problem. In 
particular in glasses and spin glasses the aging process 
displays a very rich phenomenology demanding new the- 
oretical concepts^. But already less complex — and ap- 
parently less glassy — systems, like disordered but non- 
frustrated systems^ or even pure systems^ reveal interest- 
ing and unexpected aging phenomena. One of the most 
intruiging questions in this context is whether the out- 
of-cquilibrium dynamics is essentially fully determined 
by a coarsening process (a question that even arises in 
the more complex spin glass situation—), describable by 
a growing length scale that characterizes essentially all 
out-of-equilibrium processes. In this paper we will con- 
sider a disordered system in which this question has not 
been clarified yet — and for which the answer that we 
find will reveal a non-standard scenario. 

Among glassy systems, there is a wide interest in 
disordered elastic systems, which cover a wide range 
of physical situations ranging from vortex lattices in 
superconductors^, interfaces in disordered magnets^ or 
electron glasses^ for which nonequilibrium effects are ex- 
perimentally relevant. Here, we investigate the non equi- 
librium relaxational dynamics of a solid-on-solid (SOS) 
model on a disordered substrate, defined on a two di- 
mensional square lattice and described by the following 
elastic Hamiltonian in terms of height variables hi 



SOS 



£(*< 



hi = rii + di 



(1) 



where rii are unbounded discrete variables, i.e. n% € 
{0,±1,±2, ..} and di 6 [0, 1[ are uniformly distributed 
quenched random offsets, uncorrelated from site to site. 
In the absence of disorder, i.e. di = 0, the model exhibits 
a roughening transition in the same universality class as 



the Kosterlitz-Thouless transition- , at a temperature T r 
separating a flat phase at low T from a logarithmically 
(thermally) rough one above T r . The presence of disor- 
der is known to modify significantly the nature of the 
transitio n 10 : 11 ' 12 . The so-called superroughening transi- 
tion occurs at a temperature T g = T r /2 = 2/ir. Above 
T g , where the disorder is irrelevant on large length scales, 
the surface is logarithmically rough again, although be- 
low T g the system exhibits a glassy phase where the pin- 
ning disorder induces a stronger roughness of the inter- 
face. In the continuum limit, near T g , this SOS model (JTJ 
is in the same (equilibrium) universality class as the Sine- 
Gordon model with random phase shifts, the so-called 
Cardy-Ostlund (CO) models 

H co = J d 2 x(\7ip(x)) 2 - A cos (2w(<p(x) - £(af))) (2) 

where ip(x) &] — oo,+oo[ is a continuous variable and 
£(x) £ [0, 27r[ is a uniformly distributed quenched ran- 
dom phase variable, uncorrelated from site to site, A 
being the strength of the disorder. The model @ arises 
in various contexts like the XY model in a random mag- 
netic field (without vortices) or in vortex physics where 
it describes a 2-d array of flux lines pinned by point like 
disorder—. The low temperature glassy phase (i.e. below 
T g ) of these models (^EJ) is described by a finite temper- 
ature fixed point associated with a free energy exponent 
9 = 0, which is an exact statement due to the statistical 
tilt symmetrjii^. 

Although these models have been extensively stud- 
ied, both analytically— and numericalljiiiiSiiS, these 
works have mainly focused on the equilibrium proper- 
ties. Among them the static roughness of the interface 
has been investigated thoroughly and for the dynamics 
the dynamical exponent aiiiSS. The latter was found 
to depend continuously on T and computed using the 
renormalization group (RG) up to one loop in the vicin- 
ity of T g , where the fixed point is controlled by the small 



2 



parameter r = (T g — T)/T g . Only recently, the non- 
equilibrium rclaxational dynamics (defined by a Langevin 
equation) of the Cardy-Ostlund model (J2J) was investi- 
gated analytically^ in the perturbative regime (r <C 1). 
Using the RG this study focused on the the two-times 
(t,t w ) correlation and response functions. The autocor- 
relation and local response function were found to scale 
as t/t w and characterized by asymptotically algebraic 
scaling functions with an associated decay exponent that 
depends continuously on T and was calculated perturba- 
tively up to one loop order. Finally, the associated fluctu- 
ation dissipation ratio (FDR) in the large time separation 
limit was found to be non-trivial and also T-dependent. 

In this paper we intend first to test numerically this 
analysis near T g , then to go beyond the perturbative 
regime and explore the low T dynamics where one ex- 
pects to observe a stronger signature of the logarithmic 
free energy landscape^ as suggested by the static value 
of 6 = 0. Furthermore, having determined these differ- 
ent non equilibrium dynamical properties, we propose 
to relate them to a real space analysis of the equili- 
bration process of the thermal fluctuations in the sys- 
tem. Their quantitatively precise study is possible due 
to an algorithmic*^ that allows one to compute the exact 
ground state of 

The outline of the paper is as follows. In Section II, 
we give some details of our simulations and present the 
definitions of the dynamical two times quantities we will 
focus on. In Section III, we present our numerical results 
for these quantities, and establish a comparison with the 
analytical predictions oSi (some details of this compar- 
ison are left in Appendix A). Section IV is devoted to a 
physical discussion, based an aging scenario in real space. 
Finally we draw our conclusions in the last section. 



II. SIMULATIONS AND DEFINITIONS. 

We perform a numerical study of the nonequilibrium 
relaxational dynamics of these models Q |2J) on a 2-d 
square lattice with periodic boundary conditions using 
a standard Monte Carlo algorithm. Although the SOS 
model is by definition a discrete model, the CO model 
(|2~|). which is a continuous one, needs to be discretized for 
the purpose of the simulation. We will use the discretized 
version of the gradient in with (p(x) — > tpi, i being 
the site index. The value of the displacement field ipi is 
itself discretized into 4096 intervals of width Aip between 
±4. Except when we explicitly mention it, the system is 
initially prepared in a flat initial condition (rii(t = 0) = 
or (fi(t = 0) = 0). At each time step, one site is randomly 
chosen and a move rii — > rii + 1 or ni — > rii — 1 is proposed 
with equal probability (for the CO model, the field ipi is 
incremented or decremented by an amount Aip). This 
move is then accepted or rejected according to the heat 
bath rule. Our data were obtained for a lattice of linear 
size L = 64 or L = 128, and a time unit corresponds to 
L 2 time steps. 



We will first study the connected autocorrelation func- 
tion C(t, t w ) 

C(t,t w ) = j^^{hi{t)hi{t w )) - (hi(t))(hi(t w )) , (3) 

i 

which is a two-time quantity allowing to characterize ag- 
ing properties. Then we will consider the spatial (2- 
point) connected correlation function 

= 71 E <M*)W*)> - (hidMhi+rft)) (4) 

i 

from which we measure the dynamical exponent z. In 
0J, (...) and — mean an average over the thermal noise 
and respectively over the disorder. When studying the 
CO model @ the corresponding correlation functions are 
defined by Eq. 100} with the substitution hi(t) — > tpi(t). 

These two quantities (0 QJ) are straightforwardly com- 
puted from our simulation which stores at each time step 
t the value of the height field hi(t) on each site i. Typ- 
ically, in our simulations we compute C(t,t w ) by aver- 
aging over 64 (32) different realizations of the thermal 
noise for a given configuration of the disorder and then 
averaging over 256 (128) different disorder samples for 
L = 64 (respectively L = 128). We observed that the 
main fluctuations in the computation of the correlations 
were coming from the average over the disorder. There- 
fore, we have estimated the error-bars from the sample to 
sample fluctuations of the thermal average value in (l.'il 111 . 

We are also interested in the violation of the fluctu- 
ation dissipation theorem (FDT) associated with local 
fluctuations © for which we have to consider the asso- 
ciated local linear response lZ(t,t w ) 

where fi(t w ) being an infinitesimal force applied at site 
i at time t w . The dynamical rules are then modified 
by adding a term — J^. forii to the original Hamiltonian 
Eq. ft). Numerically it is more convenient to calculate 
instead the integrated response 

p(t,t w ) = / dsK(t,s) . (6) 
Jo 

In order to isolate the diagonal component of the re- 
sponse function, we used the standard strategy22i24: we 
simulate two replicas of the system, one without an ap- 
plied force and another in which we apply a spatially ran- 
dom force to the system from time t = to time t = t w . 
This force field is of the form fi = foti, with a constant 
small amplitude /o and a quenched random modulation 
e- L = ±1 with equal probability, independently at each 
site i. The integrated response p(t,t w ) is then computed 
as 

*y4 E HEffl , f>( „, (7) 
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where (hi(t)) f i means the thermal average in the pres- 
ence of the force field /j. We have used a numerical 
value of fo — 0.3 and have checked that we were in- 
deed probing the linear response regime. Our numerical 
data for p(t,t w ) are averaged over 64 (32) independent 
thermal realizations for a given disorder configuration 
and the random fields fi for L = 64 (resp. L = 128) 
and then averaged over 512 different disorder realiza- 
tions. The error-bars are estimated in the same way 
as for the correlation functions. We point out that in- 
stead of p(t,t w ) many studies, e.g in spin glasses, focus 
on x(t]tw) — J t dsTZ(t,s). In our model in which one 
time quantities such as C(t, t) grow without bounds when 
t increases there may be a regime in which the integral 
over s in the definition of x(£, t w ) is actually dominated 
by the latest times and thus depends only very weakly 
on the waiting time t w . Therefore, in order to disentan- 
gle the off-diagonal part of the response itself the com- 
putation of p(t,t w ) JJJ), which does not suffer from the 
aforementioned peculiarity, is better suited. 

When the system is in equilibrium the dynamics is time 
translation invariant (TTI) and two-times quantities like 
C(t,t w ) or p(t,t w ) depend only on the time difference 
t — t w . Moreover, C(t,t w ) and the response lZ(t,t w ) are 
related by the Fluctuation Dissipation Ratio (FDT): 



d tw C(t,t w )=TK(t,t w ) 



(8) 



When the system is not in equilibrium, these properties 
do not hold any more and it has been proposed to gen- 
eralize the FDT to noncquilibrium situations by defining 
a Fluctuation Dissipation Ratio (FDR) X(t,t u 



a, 26. 



T 



^ d tvj C{t,t w ) 

X(t,t w ) K(t,t w ) 



(9) 



such that X(t,t w ) = 1 in equilibrium JSJ) and any devia- 
tion from unity being a signature of an out of equilibrium 
situation. In this paper, we will investigate this FDR © 
for the (noncquilibrium) rclaxational dynamics following 
a sudden quench at t = 0. Of particular interest is the 
limiting value 

X°° = lim lim X(t,t w ) , (10) 

t w — >oc t— »oo 



III. RESULTS 
A. Correlation function. 

1. Autocorrelation function. 

Figure ^ shows the decay of the connected correlation 
function C(t,t w ) for different waiting times t w and for a 
temperature T = 0.63 T g : they show a clear t w depen- 
dence. We notice that the quantity C(t w ,t w ) depends 
also on t w , before saturating to its equilibrium value for 
t w —> oo (which depends on the system size L). This 
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FIG. 1: Connected correlation function C(t,t m ) as a func- 
tion of t — t m for different waiting times t w . The inset shows 
the plot of C(t m ,t w ) — C(t,t w ) as a function of t — t w , for 
the same different waiting times, which exhibits the "quasi- 
equilibrium" regime. Here, T — 0.63 T g . 
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FIG. 2: Connected correlation function C(t,t w ) as a function 
oit/tm for different waiting times at temperature T = 0.63 T g . 
The dotted line is the result of the fit (lilt taking into account 
the data points with t/t w > 10. 



explains why one does not observe a "quasi-equilibrium" 
regime, where C(t, t w ) = C(t — t w ) when t — t w <C t w for 
the relatively small waiting times showed in Fig. ^ This 
"quasi-equilibrium" regime can however be observed if 
we plot C(t w , t w ) — C(t, t w ), as shown on the inset of Fig. 
1. 

In the aging regime, for t —t w ~ 0(t w ), these curves 
for different waiting times t w fall on a single master curve 
when we plot C(t,t w ) as a function of t/t w (Fig. |2J). In 
the large time separation regime t 3> t w these data are 
well fitted by a power law decay 



C(t,t w ) - — 



t 



-X/e 



t > U 



(11) 



Notice however that one can not exclude logarithmic cor- 
rections at low temperature where the decay exponent 
becomes very small. In Fig. [21 we plot the value of the 
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T/T g 



FIG. 3: Decay exponent X/z as a function of T/T g . The 
dashed lined indicates the exact value for T > T g . The solid 
line shows the result of the one loop RG--, given in (I I 2ft . Im- 
portantly, this curve is drawn without any fitting parameter, 
T g — 2 /it being exactly known. 



decay exponent X/z for different temperatures. In the 
high temperature phase, T > T gi where A = z = 2 one 
expects X/z = 1 independent of T (notice that the high 
temperature phase is critical and as such also displays 
aging behavio»2i2£). For T < T g the presence of disor- 
der reduces the decay exponent X/z which now depends 
continuously on temperature. In the vicinity of T g one 
observes a rather good agreement with the perturbative 
RG computation to one loop2I: 

- = l-e^T + C(r 2 ) (12) 

z 

where je = 0.577216 is the Euler constant. With the 
RG result z = 2 + 2e TB r + 0(t 2 ) this corresponds to 
A = 2 + C(r 2 ). 

Notice that the simulations near T gi T/T g > 0.8, i.e. 
in the weak disorder regime, have been performed using 
the random phase Sine-Gordon formulation @ of the 
SOS model, for which the asymptotic regime is faster 
reached for these temperatures. The inverse is of course 
true at low temperature. When it was possible, we have 
compared for a given temperature the asymptotic prop- 
erties of C(t,t w ) using the SOS model and the CO 
model (J5J. We show the result of this comparison for 
T = 0.63 T g in Fig. H 

One observes that both formulations are in good agree- 
ment concerning the t/t w scaling form and are in a rea- 
sonable agreement concerning the value of the exponent 
X/z, thus confirming the universality of this property. 
However, the amplitude itself does not seem to be uni- 
versal. 

At lower temperature the perturbative calculation fails 
to predict the correct behavior of X/z: in Fig. |3we ob- 
serve a change in its T-dependence below T w 0.8 T g . In 
this regime one obtains a good fit of the decay exponent 



FIG. 4: Connected correlation function C(t,t w ) obtained with 
the SOS model (filled symbols) and with the CO hamiltonian 
(open symbols) as a function of t/t w for different t w . Here 
T = 0.63 T g . 

by 

~~A X/Z T , A x/Z =0.85 ±0.04 (13) 

If one naively assumes that the one loop RG calculation 
A = 2 is still valid at low temperature, this would already 
indicate a 1/T behavior of the dynamical exponent z. We 
will come later to this point where we explicitly compute 
this exponent z. Indeed, this scaling form (|11|) can be 
written as 

^-(z^r ' (14) 

thus defining a length scale C(t) which can be further 
analyzed by measuring how the spatial correlations are 
growing in the system (see the next paragraph). The 
functional shape of C(t,t w ) that we determined suggests 
that its T-dependence is mainly contained in the decay 
exponent within the the aging regime where (t — t w ) ~ 
0(t w ). It is remarkable that its most prominent feature, 
the t/t w scaling and the asymptotically algebraic scal- 
ing form with a T-dependent decay exponent, is already 
captured by the one loop RG calculation oSi. By con- 
trast, one observes that the "quasi-equilibrium" regime 
(t — t w ) <C t w shows a much stronger T-dependence. 
At low temperature T < T g /2 the autocorrelation func- 
tion C(t,t w ) displays an inflection point at small time 
difference t — t w . In Fig. [31 where C(t,t w ) as a func- 
tion of t — t w is shown in a linear-log plot for differ- 
ent large waiting times t w , one observes a qualitative 
change of behavior which could suggest a finite limiting 
value lim t _ (00 limt^^oo C(i, tu,). However, on the time 
scales explored here, we have not identified a clear sig- 
nature of such a behavior. Nevertheless, this point de- 
serves further investigations of the equilibrium properties 
at low temperature, where some discrepancies between 
numericsi&iS and analytical predictions^ were already 
found. 
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FIG. 5: Autocorrelation function C(t,t w ) as a function of 
t — t w for different large waiting times t w , at (very) low tem- 
perature, T — 0.15 T g . For short t — t w , this quantity shows 
an inflection point. The inset shows the same quantity for 
T = 0.63 T g , which exhibits a qualitatively different behavior 
for t — tw <U. 




FIG. 6: Spatial connected correlation function C(r,t) as a 
function of r for different times t. Here T — 0.47 T„. 



2. 2-point correlation function. 

In Fig. [5]we show the 2— point correlation function (JIJ 
for a temperature T = 0.47 T g (and L = 64) for different 
times t. 

As t grows spatial correlations develop in the system. 
More precisely, as shown in Fig. [7| C(r, t) scales as 



C(r, t)=T 



C(t) - 1 1 ' 7 - 



(15) 



The value of z that gives the best data collapse leads to 
our first estimate of the dynamical exponent. The log- 
arithmic behavior for r <C C{t), C(r,t) ~ lnX(t)/r is 
in agreement with the constraint imposed by the statis- 
tical tilt symmetry (STS)i£ which fixes the equilibrium 



FIG. 7: Spatial connected correlation function C(r,t) as a 
function of r/t 1/lz with 1/z = 0.17 ± 0.01 for different times 
t. Here T = 0.47 T„ 



behavior of the connected 2-point correlation function to 

2 T 



lim C(r, t) 



lnr 



(16) 



which is identical with the pure (i.e. disorder-free be- 
havior). We also checked that the amplitude of the log- 
arithmic behavior of C(r,t) for r/C(t) <C 1 is in good 
agreement (within a few percents) with Eq. (|16|) . 



3. Dynamical exponent. 

Another way to estimate the dynamical exponent is 
to determine the time dependent length scale C(t) itself. 
For that purpose, and given the scaling form previously 
computed (| 1 51) . we estimate C(t) via a the space integral 
of the spatial correlations^: 



L/2 



L/2 



drC(r,t)= drT{r/£(t))~C{t) duT(u), 
Jo Jo 

(17) 

where we assumed in the last step that L/C(t) <C 1 
(which is indeed the case on the time scales considered 
here) and that C(r,t) decays sufficiently fast at large r 
(we checked that it actually decays exponentially). No- 
tice also that the sum in (|17l) is bounded to L/2 due 
to periodic boundary conditions. In Fig. [S] we showed 
the value of C(t) computed with (|17|l for different tem- 
peratures. One obtains a rather good fit of these curves 
(Fig. Q by a power law C(t) ~ t 1 / 2 ^, thus obtain- 
ing a value of the T dependent dynamical exponent in 
good agreement with the value obtained by collapsing 
the different curves in Fig. [7| One notices also that C(t) 
approaches an algebraic growth after a pre-asymptotic 
regime which increases with decreasing temperature. Fig. 
[5] shows our estimate for l/z(T) as a function of T. As 
expected, the dynamical exponent is a decreasing func- 
tion of the temperature. One expects that z — 2 for 
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FIG. 8: Growing length scale £(t) computed from 1171 for 
different temperatures. The solid lines are guides to the eyes. 



FIG. 10: Integrated response function p(t,t w ) as a function 
oit — tw for different waiting times t w . Here T = 0.47 T 9 . 
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FIG. 9: l/z(T) as a function of T/T a . The dashed line which 
shows the result of the one loop R G 11,20 is drawn without any 
fitting parameter. 



T > T g and that it becomes T-dependent below T g with 
z = 2 + 2eJ E r + 0(t 2 ) as predicted by a one loop RG 
calculation 11,20 . At high temperature T > T g and in the 
vicinity of T~ , it is numerically rather difficult to extract 
a reliable estimate for the dynamical exponent from 115|) 
or (|17|) due to finite size effects. Therefore we restrict 
ourselves here to lower temperatures T < 0.8T g . For 
temperature T > 0.7 T g , the value of z is still in rea- 
sonable agreement with the RG prediction. Around the 
value T* ~ 0.63 T g , where z ~ 4, the curve l/z(T) shows 
an inflection point, below which 1/z decreases linearly 
with T. In this regime, z(T) is well fitted by 



z(T) 



forT < T* 



(18) 



which, given (|13fl . shows also that A ~ 2 is still a good 
estimate a low T. This behavior z oc l/T is compati- 
ble with an activated dynamics over logarithmic barri- 
ers, i.e. an Arrhenius type behavior t typ ~ e B ^t yP / T w ith 
Bl ~ logLtyp- Assuming that the largest barriers, 
which dominate the low temperature dynamics, encoun- 



tered in this non equilibrium relaxation process have the 
same scaling that the equilibrium ones, this logarithmic 
behavior is also consistent with a free energy exponent 
6 = C 32 . Interestingly, this change of behavior of z at 
a value of z c = 4 above which z oc l/T (|18fl is reminis- 
cent of the related case of a particle in a one dimensional 
disordered potential with logarithmic correlations where 
such a behavior was obtained analytically 22 . It should 
be mentioned that a dynamical exponent that varies like 
1 /T has also be found in other disordered systems like in 
spin glasses^2i2& and in random ferromagnets^i. Finally, 
although 118fl suggests the existence of a well defined typ- 
ical relaxation time, one expects the full distribution of 
the barrier heights to be very broad^ 3 -, and needs proba- 
bly further work to be investigated. 



B. Integrated response function. 

In this section, we focus on the integrated response JTJl . 
In Fig. ijTU|l we show a plot of p(t,t w ) as a function of 
the time difference t — t w for different waiting times t w . 
Here too, one observes a clear waiting time dependence. 

These curves for different waiting times t w fall on a 
single master curve if one plots them as a function of 
t/t w , as shown in Fig. ^2 As suggested on this log- 
log plot (Fig. [TTjl . p(t,t w ) takes the following power law 
decay 



p(t,t w ) 



tu 



-A/2 



t > t v 



(19) 



Notice that the decay exponent, within the accuracy of 
the data presented here, is the same as the one of the cor- 
responding autocorrelation function C(t,t w ) l|ll|l. This 
t/t w scaling form, together with the relation between the 
decay exponent of p(t, t w ) and C(t, t w ) are also fully com- 
patible with previous one loop RG calculations. As we 
will see, this has important implications for the FDR as 
discussed in the next paragraph. 
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FIG. 11: Integrated response function p(t,t w ) as a function 
of t/t w for different waiting times t m at T = 0.47T ff . 



C. Fluctuation dissipation ratio. 

In order to characterize the deviation from the equilib- 
rium, we compute in this section the FDR X(t,t w ) 10- 
For T > T g the disorder is irrelevant and the FDR is ex- 
pected to be identical to the FDR of the pure case, which 
we therefore consider first: in the pure model the aucor- 
relation and the response function can be computed an- 
alytically. In the out of equilibrium regime t w < t <C L 2 
(remembering that z = 2 for the pure case) one ha& 

1 1 

>?>pxae(t, t w ) — . vnT 7 j t > t w 
1 g\ 2lr ) t ~t w 

<W*'*-) = iT^ ta (j£%) (20) 

Using these expressions (I2UI1 together with JjJJ, one ob- 
tains that X(t,t w ) = X(C(t,t w )) which allows to write 
the relation defining the FDR 10 in an integrated form 
using the definition of p(t,t w ) © 

,*„ 

T Ppure {t : t w ) / dsX purc (^purc {t : s) )d s Cp Ure (t, s) 
JO 

— Xp Uie (C(t,t w )) — X puxe (C puxe (t, 0)) (21) 

with d u X purc (u) = X pure (u). C purc (t, 0) is expected to be 
small one can extract X purc (C(t,t w )) from the slope of 
the curve Tp purc (t, t w ) versus C puie (t, t w ) in a parametric 
plot, provided t w is sufficiently large such that the curves 
for different t w collapse. In Fig. 1)12(1 this parametric plot 
Tppure versus C puro is shown. For large values of C pure one 
expects to recover the FDT, and a slope of value unity. 
On the other hand, as C purc decreases all these curves 
converge to a same master curve X pure (C), which, using 
(|2"U1) can be exactly computed for the pure model 

c 

«7)=7ln^ , (22) 



As one can see in Fig. our numerical results are 
in good agreement with the exact calculation. An im- 
portant point is that the slope at the origin gives the 
asymptotic value of the FDR X?^ e , Eq. (fTUf) such that 

Tppurcit, tw*) ^ -^rjure ' ^ purey ^s tw) when C pure (t , t w ) ► 0. 

As is obvious from Eq. ((22(1 for the pure model one has 
^pure = 1/2; the random walk valued, independent of the 
temperature. 

For a finite size system, one expects to recover the 
equilibrium dynamical regime for large but finite waiting 
times t w 3> tEQ and in particular the restoration of the 
FDT © reflected by X(t,t w ) = 1. Therefore, as pre- 
dicted by the analytical solution, the parametric curve 
of Tp vs. C will progressively move to the right with 
increasing t w converging in equilibrium (t w —* oo) to a 
straight line passing through the origin. 

We now turn to the case T < T g when the disorder is 
relevant. Given the t/t w scaling forms we have obtained 
for C(t, t w ) (Fig. |2J) and for p(t, t w ) (Fig. ITTjl one expects 
also in the disordered case to have p(t, t w ) = X(C(t, t w )). 
Indeed, as shown in Fig. 1131 the parametric plot Tp vs. 
C for different t w is qualitatively similar to the curve 
obtained for the pure case. In particular, the property 
p{t,t w ) = X(C(t,t w )), together with Eq. yields, in 
the nonequilibrium regime 

Moreover, our data (Fig. 113(1 are consistent with a finite 
limiting value (as defined in Ea. l(lU(l ) X°° > also in 
presence of disorder although the asymptotic value of this 
quantity is very difficult to estimate numerically. This 
fact is qualitatively in agreement with RG predictions. 
In contrast to the pure model, and according to—, this 
value X°° depends continuously on T as 

X°° = - + 0{t 2 ) (24) 

z 

close to Although a precise comparison with this 

RG predictions in the vicinity of T g , where the deviations 
from the pure case are expected to be small, is difficult at 
this stage (requiring a study on longer time scales) one 
can see in Fig. ^| and Fig. that our data are still in 
a reasonable agreement with the one loop relation ((24(1 . 

IV. COARSENING OR GROWING 
FLUCTUATIONS? 

The behavior we obtained for the 2-point correlation 
function C(r, t) allowed to identify a growing length scale 
C(t) on which the system gets equilibrated. To go fur- 
ther, one would like to relate this length scale £{t) to 
the size of spatially correlated structures like domains or 
droplets. We first explored the idea that at low temper- 
ature, the nonequilibrium dynamics could be understood 
as a coarsening process reflected in a spatially grow- 
ing correlation with the ground state (GS). Interestingly, 
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FIG. 12: Parametric plot of the integrated response function 
Tp(t,t w ) as a function of C(t,t w ) for different waiting times 
tw and T — 1.1 T g . The solid line is the result for the pure 
case as given by Eq. [122ti and does not contain any fitting 
parameter. The dashed line shows the slope corresponding to 
the non-violated FDT. 




FIG. 13: Parametric plot of the integrated response function 
Tp(t,t w ) as a function of C(t,t w ) for different waiting times 
t w . Here T = 0.47 T g . The solid line corresponds to a value 
of X°° — 1/z 12411 although the dashed one corresponds to 
Xoo = 1/2, thus showing a clear deviation from the pure 
case. The dotted line shows the slope corresponding to FDT. 



computing the GS of the SOS model on a disordered sub- 
strate Q is a minimum cost flow problem for which exits 
a polynomial algorithm and can therefore be computed 
exactljsiSiiS. After determining one GS (note that the 
GS, which is computed with free boundary conditions, is 
infinitely degenerated since a global shift of all heights 
by an arbitrary integer is again a GS), we define for each 
time t the height difference mj(i) = rii(t) — rii(0) and 
identify the connected clusters (domains) of sites with 
identical rrn(t) using a depth- first search algorithm. No- 
tice that for the comparison with the Ground State, the 
Monte Carlo simulations are performed here using free 
boundary conditions. 

In Fig. 1151 we show snapshots of these domains for 



FIG. 14: Comparison between X°° , open symbols and 1/z, 
filled symbols. The value of 1/z for T = 1.1 T g shown here is 
the exact one. 



T > T g on the left panel and T < T g on the right one. 

Starting from a random initial configuration one can 
for T < T g very quickly (t < 100) identify large domains 
that evolve only very slowly at later times. On the other 
hand for T > T g the configurations decorrelated very 
quickly in time. To make this analysis more quantita- 
tive, we determined the cluster size distribution P t ^(S,t) 
for one realization of the disorder (and for different real- 
izations of the thermal noise). 

As shown on Fig. [23 Pth(S, t) starts to develop a peak 
at a rather large value S*(t) on the earlier stage of the 
dynamics (this peaks also develops if we start with a ran- 
dom initial configuration). It turns out that S*(t) is the 
size of the largest connected flat cluster of the ground 
state configuration n® = C st . On the time scales pre- 
sented here, as time t is growing, this peak remains stable 
S* (t) ~ C st , implying that the system is not coarsening. 
At later times, as suggested by simulations on smaller 
systems, this peak progressively disappears and the dis- 
tribution becomes very flat. We also checked that the 
mean size of these connected clusters is not directly re- 
lated to C(t). 

One has however to keep in mind that we are comput- 
ing the connected correlation functions, i.e. we measure 
the thermal fluctuations of the height profile around its 
mean (typical) value (hi(t)). Therefore, we believe that 
these connected correlations are instead related to the 
broadening of this "stable" peak (Fig. I16fl . i.e. the fluc- 
tuations around this typical state at time t. The slow evo- 
lution of the typical configuration, compared to the one 
of thermal fluctuations around it, is corro borated by the 
one loop calculationSiiM which shows that (hi(t)) (hi(t w )) 
decays as 



(hi(t)){hi{t w )) ~r ( i- 



-1/2 



+ G(r 2 ), 



(25) 



i.e. much slower than the connected one l|Sl ll2|) . 

To characterize more precisely the fluctuations of this 
cluster, we have followed the following protocol: after a 
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FIG. 15: Snapshot of the height field relative to the ground state rrii(t) = ni(t) — n® for T > T g on the left panel and T = 0.47 T g 
on the right panel. Different colors correspond to different values of rrii(t) : mi(t) = —2 in green, rm(t) = — 1 in white, rrn(t) = 
in black and mi(t) = +1 in blue and so on. Note that for T > T g the configuration at t — 10° is already decorrelated from the 
one at t = 10 3 , whereas forT < T g large domains in white and black persist and change only slowly in time. 



time ti ~ 100 we store the configuration of the largest 
connected cluster. Then, for each time t, we compute 
the distribution ^droplet °f the size of the connected 
clusters that were part of this cluster at time U but not 
at time t (the subscript "flat" refers to the flat initial 
condition). In Fig. [T7| we show a plot of ffroplet(^) *) 



for a temperature T = 0.47 T g , for different times t. 

It decays as a power law for small sizes S, and this 
power law behavior extends to larger and larger values 
of S as t is growing. Although these data give already 
some interesting insight on how the thermal fluctuations 
equilibrate in the system, it turns out to be very hard to 
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FIG. 16: Size distribution Pth(S, i) (see definition in the text) 
for different times t. Here T = 0.47T o . 



FIG. 18: Distribution of the size of the clusters ^dropietCS*! i) 
as a function of S and for different times t. Here, the initial 
condition is the Ground State and T = 0.3 T g . 
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FIG. 17: Size distribution Pdropiet( S > *) (symbols) for differ- 
ent times t. The solid lines represent fdropiet *) ( see the 
definition in the text) at the same corresponding times. Here 
T = 0.47 T„. 



FIG. 19: S a P^ plct (S,t) with a = 1.9 ± 0.1 as a function of 
S/t 2/z with 2/z = 0.26 ± 0.03. Here the initial condition is 
the Ground State and T = 0.3 T„. 



obtain good statistics for larger values of S in this way. 
In order to perform a more precise quantitative anal- 
ysis of this distribution we identify alternatively these 
"droplets" by initializing the system in the ground state 
itself Rj(t = 0) = n®. At low temperature, and on the 
time scales explored here, one expects that the ground 
state represents a good approximation of a typical con- 



figuration, i.e. (ni(t)} 



Again we compute the 



distribution P '^piet (S, t) of the sizes of the connected 
clusters with a common value of rrii(t) ^ 0. As shown 
in Fig. El-PdJopietC 1 ^'*) determined in this way coincides 
very well with Pfroplet f )- Moreover, the calculation of 
^droplet i s mu ch easier and allows for a more precise 
analysis. 

In Fig. we show a plot of -P^opiot extending to 
larger times, for a temperature T — 0.3 T g . It turns out, 
as shown in Fig. El that -P^opict obeys the scaling 



form 



7 GS 



s 



, a = 1.9 ±0.1 
(26) 

where a is independent of T within the accuracy of our 
data and C{t) ~ i 1 / 2 . The value of z in l|2()l) is in good 
agreement with the one extracted from the 2-point cor- 
relation function C(r,t) = F(r/C(t)) (|15fl . Furthermore, 
considering that each "droplet" of size S > r 2 gives a 
contribution to C(r,t) proportional to S, one obtains, 
given the distribution (12(:i[l with a = 2 



C(r,t) oc 
oc In C(t)/r 



dSSP^ plct (S, t) 
, C(t)/r < 1 



(27) 



which is consistent with the behavior we obtained in Fig. 
and Eq. (Ify 

This scaling form (|26|l thus establishes a relation be- 
tween £(t) and the typical size of compact excitation 
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around a "typical" configuration, evolving slowlier in 
time. 



V. CONCLUSION. 

In conclusion, we have performed a rather detailed 
analysis of the nonequilibrium relaxational dynamics of 
the SOS model on a disordered substrate (JTJ, and of the 
related Cardy Ostlund model Close to the super- 
roughening temperature T g our results for the autocorre- 
lations, spatial correlations and response function as well 
as for the fluctuation dissipation ratio (FDR) agree well 
with the prediction of a recent one loop RG calculation 21 , 
whereas deep in the glassy low temperature phase sub- 
stantial deviations occur. 

The aging features obtained perturbatively, charac- 
terized by a t/t w scaling of local correlation and re- 
sponse functions with a temperature dependent decay 
exponent, carries over into the low temperature regime, 
including a non-trivial temperature dependent fluctua- 
tion dissipation ratio Xoo associated to these correlation 
and response functions. The change in the low temper- 
ature behavior of these quantities compared with the 
RG predictions turns out to be contained in a change 
of the functional temperature dependence of the dynam- 
ical exponent z(T), which relates the age t of the system 
with a length scale C{t): z(T) changes from a linear In- 
dependence close to T g to a 1/T-behavior far away from 
T g . This is, to our knowledge, the first clear indication 
of an activated dynamics over logarithmic barriers in this 
marginal glass phase {i.e. 9 — 0). Given the strong sim- 
ilarity of the behavior of z with the one found for the 
related model of a particle in a 1-d disordered poten- 
tial with logarithmic correlations^, an open question re- 
mains whether this dynamical crossover admits a static 
counterpart as found in that model 2 ^. 

The growing length scale C(t), increasing algebraically 
with the age of the system, turned out to be connected 
to the typical size of the fluctuations around metastable 
configurations with long life time in which the system 
gets trapped immediately after a quench into the low 
temperature phase. In contrast to a standard coarsen- 
ing process, where the growing length scale represents 
the typical size of domains (which are identified as spa- 
tial regions strongly correlated with one of the ground 
states of the systems) we encounter here a scenario in 
which already soon after a temperature quench theses 
domains are actually very large, but do not grow fur- 



ther and are destroyed by fluctuations of increasing spa- 
tial extent. Moreover, these fluctuations itself can be 
again identified as connected patches of ground state, or 
droplets. 

The emerging picture for the aging dynamics below the 
superroughening transition within the glassy low temper- 
ature phase thus differs from various well established ag- 
ing sceanrios in glasses, spin glasses and other disordered 
systems: As pointed out above the approach to equilib- 
rium is not a coarsening process as it occurs in other dis- 
ordered systems like the random ferromagnet£. It also 
differs from the aging process encountered in finite di- 
mensional spin glasses which also display coarsening 2 *^ 
with doamins that can straighforwardly be identified due 
to the existence of the Edwards- Anderson order param- 
eter. On the other hand the aging scenario revealed for 
this system appears to be far from being as complex as 
in mean-field spin glasses^: It is more reminiscent of the 
dynamics of a random walk in a one-dimensional energy 
landscape, the Sinai-model, in which the walker displace- 
ment also increases only logarithmically with time due to 
the existence of deep traps with exponentially long trap- 
ping tiniest. 

With regard to our observation that these traps in the 
disordered SOS model can be identified with configura- 
tions roughly made of large patches of the ground state 
it is tempting to describe the aging process here as a dif- 
fusion in a coarse grained configuration space consisting 
of height profiles composed like a jigsaw puzzle of ground 
state domains of optimzed shape (most probably flat 
pieces of constant height with energy minimizing bound- 
aries). The escape from a deep energy minima proceed, 
according to what our numerical analysis indicates, via 
the thermal activation of larger and larger patches, each 
intermediate configuration again being metastable with 
some finite survival time. This process is reminiscent of 
the energy- well- within-energy- well picture proposed in«& 
and in our view further studies would be worthwhile to 
develop this analogy in more detail. 
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inal and the ones in real space computed numerically 
in the present paper. We give here the details for the 
connected autocorrelation function C(t,t w ) the ex- 
tension to the integrated response p(t, t w ) (|7J) being then 
straighforward. In2i, the analytical predictions focused 
on the following connected correlation function 



C q {t,t w ) = (h q (t)h- q (t w )) - (h q (t))(h- q (t w )) (Al) 

where h q {t) is the Fourier transform, w.r.t. space vari- 
able, of the field hi(t) Using RG along the line of 
fixed points near T g , this correlation function 1|A1|) was 
computed up to order C(r 2 ). It takes the following form 



T ( t 



Ci(t,t w ) = ^( — ) F c {q z {t~t w ),t/t w ){k2) 

q 2 \t w J 

c = e IE T + O(T 2 ) 

where je is the Euler constant, given in the text, and 
with the asymptotic behavior in the large time separation 
limit 



F c {v,u) = \-0{u ) 



(A3) 



The connected autocorrelation function C(t,t w ) @ we 
compute here is related to <|A1|1 through 



d 2 q 



Sc r .12 



(t , t w ^j 



C(t, t w ^ - 



(A4) 



Perfoming the change of variable v — q z (t — t w ), (|A4|) 
becomes 

T ( t \ ^ c f°° dv 
C(t,t w ) = — ( — ) / —F c (v,t/t w ) (A5) 

27TZ \t w J J V 

where we have taken the IR (resp. the UV) cutoff to 
(resp. to oo) and checked the convergence of the integral 
over v. Using the asymptotic behavior l)A3|l one obtains 
(the remaining integral over v being well defined) in the 
large time separation limit t 3> t w 



T ( t 

c(t,t w )~ — (- 

Zirz \ ti„ 



9c-l ^oo 



/ -F Coo (v) (A6) 
Jo v 



APPENDIX A: COMPARISON WITH RG 
CALCULATIONS NEAR T g 

In this appendix we establish the connexion between 
the quantities (in Fourier space) computed analytically 



which, given the value of 9c l|A2|l . leads to the following 
one loop result for the decay exponent X/z <|ll|l : 



X/z= 1 - e 7E r + 0(t 2 ) 
given in the text in Ea. (|12fl . 



(A7) 



